
/*------------------------------------------------------------------------*/
/*  File: cluster.c                                                       */
/*  Version 1.0                                                           */
/*------------------------------------------------------------------------*/
/*  Written for Software Project Class.                                   */
/*  Based uopn qpex1.c                                                    */
/*  Permission is expressly granted to use this example in the            */
/*  course of developing applications that use ILOG products.             */
/*------------------------------------------------------------------------*/

/* cluster.c - Entering and optimizing k-clustering integer linear programming problem */

#include "cluster.h"

int kMax;
int nVertices;
int nEdges;

/* This routine initializes the cplex enviorement, sets screen as an output for cplex errors and notifications,
 and sets parameters for cplex. It calls for a mixed integer program solution and frees the environment.
 To Do:
 Declare the parameters for the problem and fill them accordingly. After creating the program thus copy it into cplex.
 Define any integer or binary variables as needed, and change their type before the call to CPXmipopt to solve problem.
 Read solution (both objective function value and variables assignment). Add out parameter to the calling function to
 return the objective function max solution.
 */
int k_cluster(int inputKMax, double **xx, double *objvall)
{
    
    /* Declare pointers for the variables and arrays that will contain
     the data which define the LP problem. */
    
    char *probname = NULL;
    int numcols;
    int numrows;
    int numvars;
    int objsen;
    double *obj = NULL;
    double *rhs = NULL;
    char *sense = NULL;
    int *matbeg = NULL;
    int *matcnt = NULL;
    int *matind = NULL;
    double *matval = NULL;
    double *lb = NULL;
    double *ub = NULL;
    char *ctype = NULL;
    int *indices = NULL;
    
    int solstat;
    double objval;
    double *x;
    
    /* Declare and allocate space for the variables and arrays where we
     will store the optimization results including the status, objective
     value and variable values. */
    
    CPXENVptr env              = NULL;
    CPXLPptr  lp               = NULL;
    int       status=0;
    char buffer[512];
    
    
    nVertices = getNumVertices();
    
    nEdges = getNumEdges();
    
    kMax = inputKMax;
    
    
    
    /* Initialize the CPLEX environment */
    env = CPXopenCPLEX (&status);
    if ( env == NULL ) {
        char  errmsg[1024];
        fprintf (stderr, "Error: Could not open CPLEX environment.\n");
        CPXgeterrorstring (env, status, errmsg);
        fprintf (stderr, "%s", errmsg);
        goto TERMINATE;
    }
    
    /* Turn on output to the screen */
    status = CPXsetintparam (env, CPX_PARAM_SCRIND, CPX_ON); 
    if ( status ) {
        fprintf (stderr, "Error: Failure to turn on screen indicator, error %d.\n", status);
        goto TERMINATE;
    }
    
    /* Fill in the data for the problem.  */
    status = setProblemData (&probname, &numcols, &numrows, &numvars, &objsen, &obj,
                             &rhs, &sense, &matbeg, &matcnt, &matind, &matval, 
                             &lb, &ub, &ctype, &indices);
    if ( status ) {
        fprintf (stderr, "Failed to build problem data arrays.\n");
        goto TERMINATE;
    }
    
    
    /* Create the problem. */
    lp = CPXcreateprob (env, &status, probname);
    if ( lp == NULL ) {
        fprintf (stderr, "Error: Failed to create problem.\n");
        goto TERMINATE;
    }
    
    /* Now copy the problem data into the lp */
    status = CPXcopylp (env, lp, numcols, numrows, objsen, obj, rhs,
                        sense, matbeg, matcnt, matind, matval,
                        lb, ub, NULL);
    if ( status ) {
        fprintf (stderr, "Failed to copy problem data.\n");
        goto TERMINATE;
    }
    
    /* Now copy the ctype array */
    status = CPXchgctype(env, lp, numcols, indices, ctype);
    if ( status ) {
        fprintf (stderr, "Failed to change ctype\n");
        goto TERMINATE;
    }
    
    
    /* Optimize the problem. */
    status = CPXmipopt (env, lp);
    if ( status ) {
        fprintf (stderr, "Error: Failed to optimize problem.\n");
        goto TERMINATE;
    }
    
    /* Finally, write a copy of the problem to a file. 
    status = CPXwriteprob (env, lp, probname, NULL);
    if ( status ) {
        fprintf (stderr, "Error: Failed to write LP to disk.\n");
        goto TERMINATE;
    }*/
    
    
    solstat = CPXgetstat (env, lp);
    x = (double*)calloc(numcols, sizeof(double));
    
    /* Obtain solution */
    status = CPXsolution(env, lp, &solstat, &objval, x, NULL, NULL, NULL);
    if ( status ) {
        fprintf (stderr,"No Solution.  Exiting...\n");
        goto TERMINATE;
    }
    
    /*
     printf ("\nSolution value  = %f\n\n", objval);
     printf ("\nSolution status = %d\n", solstat);
     
     CPXgetstatstring (env, solstat, buffer);
     printf ("\nSolution status is %s\n", buffer);
     */
    
    
    *xx = x;
    
    *objvall = objval;
    
TERMINATE:
    
    /* Free up the problem as allocated by CPXcreateprob, if necessary */
    if ( lp != NULL ) {
        status = CPXfreeprob (env, &lp);
        if ( status ) {
            fprintf (stderr, "Error: CPXfreeprob failed, error code %d.\n", status);
        }
    }
    
    /* Free up the CPLEX environment, if necessary */
    if ( env != NULL ) {
        status = CPXcloseCPLEX (&env);
        
        /* Note that CPXcloseCPLEX produces no output,
         so the only way to see the cause of the error is to use
         CPXgeterrorstring.  For other CPLEX routines, the errors will
         be seen if the CPX_PARAM_SCRIND indicator is set to CPX_ON. */
        
        if ( status ) {
            char  errmsg[1024];
            fprintf (stderr, "Could not close CPLEX environment.\n");
            CPXgeterrorstring (env, status, errmsg);
            fprintf (stderr, "%s", errmsg);
        }
    }
    
    /* Free up the problem data arrays, if necessary. */
    free_and_null ((char **) &probname);
    free_and_null ((char **) &obj);
    free_and_null ((char **) &rhs);
    free_and_null ((char **) &sense);
    free_and_null ((char **) &matbeg);
    free_and_null ((char **) &matcnt);
    free_and_null ((char **) &matind);
    free_and_null ((char **) &matval);
    free_and_null ((char **) &lb);
    free_and_null ((char **) &ub);
    free_and_null ((char **) &ctype);
    free_and_null ((char **) &indices);
    
    
    return (status);
}  


/* This simple routine frees up the pointer *ptr, and sets *ptr to NULL */
void free_and_null (char **ptr)
{
    if ( *ptr != NULL ) {
        free (*ptr);
        *ptr = NULL;
    }
} 


/* MY IMP */

int setProblemData (char **probname_p, int *numcols_p, int *numrows_p, int *numvars_p,
                    int *objsen_p, double **obj_p, double **rhs_p,
                    char **sense_p, int **matbeg_p, int **matcnt_p, 
                    int **matind_p, double **matval_p,
                    double **lb_p, double **ub_p, char **ctype_p, int **indices_p)
{
    
    char     *zprobName = NULL;     /* Problem name <= 32 characters */        
    double   *zobj      = NULL;
    double   *zrhs      = NULL;
    char     *zsns      = NULL;
    int      *zmatBeg   = NULL;
    int      *zmatCnt   = NULL;
    int      *zmatInd   = NULL;
    double   *zmatVal   = NULL;
    double   *zlb       = NULL;
    double   *zub       = NULL;
    char     *zctype    = NULL;
    int      *zindices  = NULL;
   
    
    double   *protoObj    = NULL;
    char     *protoSns    = NULL;
    double   *protoRhs    = NULL;
    int      *protoMatBeg = NULL;
    int      *protoMatCnt = NULL;
    int      *protoMatInd = NULL;
    double   *protoMatVal = NULL;
    
    int nzConstI=2, nzConstII=1, status=0;
    int k, i, x=0;
    
    int nzVarsPerCluster = nEdges * (nzConstI+nzConstII);
    int nxVarsPerCluster = nEdges * (nzConstI+nzConstII*2) + 2*nVertices;
    int nColsPerCluster  = nVertices + nEdges;
    int nRowsPerCluster  = nzVarsPerCluster + nVertices + 1;
    int nVarsPerCluster  = nzVarsPerCluster + nxVarsPerCluster;
    
    int nCols = (nColsPerCluster)*kMax;
    int nRows = (nzVarsPerCluster+1)*kMax + nVertices;
    int nVars = (nVarsPerCluster)*kMax;
    
    
    zobj      = (double*)calloc(nCols, sizeof(double));
    zsns      = (char*)calloc(nRows, sizeof(char));
    zrhs      = (double*)calloc(nRows, sizeof(double));
    zmatBeg   = (int*)calloc(nCols, sizeof(int));
    zmatCnt   = (int*)calloc(nCols, sizeof(int));
    zmatInd   = (int*)calloc(nVars, sizeof(int));
    zmatVal   = (double*)calloc(nVars, sizeof(double));
    
    zprobName = (char *)malloc(32 * sizeof(char)); 
    zlb       = (double*)calloc(nCols, sizeof(double));
    zub       = (double*)calloc(nCols, sizeof(double));
    zctype    = (char*)calloc(nCols, sizeof(char));
    zindices  = (int*)calloc(nCols, sizeof(int));
    
    if (zprobName == NULL || zobj    == NULL ||
        zrhs      == NULL || zsns  == NULL   ||
        zmatBeg   == NULL || zmatCnt == NULL ||
        zmatInd   == NULL || zmatVal == NULL ||
        zlb       == NULL || zub     == NULL ||
        zctype    == NULL || zindices == NULL)  {
        
        status = 1;
        goto TERMINATE;
    }
    
    
    sprintf(zprobName, "k[%d]_cluster.lp", kMax);
    
    status = createProtoType(&protoObj, &protoMatBeg, &protoMatCnt, &protoMatInd, &protoMatVal, &protoSns, &protoRhs);
    if ( status ) {
        fprintf (stderr,"Cant create protoType...\n");
        goto TERMINATE;
    }
    
    /* duplicating protoType for each k */
    for(k=0; k<kMax; k++){
        
        /* zobj, zmatBeg, zmatCnt */
        for(i=0; i<nColsPerCluster; i++){
            
            zobj[k*(nColsPerCluster) + i] = protoObj[i];
            
            zmatBeg[k*(nColsPerCluster) + i] = protoMatBeg[i] + k*(nVarsPerCluster);
            
            zmatCnt[k*(nColsPerCluster) + i] = protoMatCnt[i];
        }
        
        /* matInd, matVal */
        for(i=0; i<nVarsPerCluster; i++){
            
            if (protoMatInd[i] < 0){
                
                if(protoMatInd[i] == -3) {
                    zmatInd[k*(nVarsPerCluster) + i] = nzVarsPerCluster*kMax + x%nVertices;
                    x++;
                }
                
                if(protoMatInd[i] == -4) {
                    zmatInd[k*(nVarsPerCluster) + i] = nzVarsPerCluster*kMax + nVertices + k;
                }               
            }
            else{
                
                zmatInd[k*(nVarsPerCluster) + i] = protoMatInd[i] + k*nzVarsPerCluster;
            }
            
            zmatVal[k*(nVarsPerCluster) + i] = protoMatVal[i];
        }
        
        /* sns, rhs */
        for(i=0; i<nRowsPerCluster; i++){
            
            if(protoSns[i] < 0){
                
                zsns[i-nzVarsPerCluster + nzVarsPerCluster*kMax] = 0-protoSns[i];
                zrhs[i-nzVarsPerCluster + nzVarsPerCluster*kMax] = protoRhs[i];
                
                if(i == nRowsPerCluster-1){
                    zsns[i-nzVarsPerCluster + nzVarsPerCluster*kMax + k] = 0-protoSns[i];
                    zrhs[i-nzVarsPerCluster + nzVarsPerCluster*kMax + k] = protoRhs[i];
                }
            }
            else{
                
                zsns[k*nzVarsPerCluster + i] = protoSns[i];
                zrhs[k*nzVarsPerCluster + i] = protoRhs[i];
            }
        }
    }
    
    for(i=0; i<nCols; i++){
        
        zlb[i] = 0;
        zub[i] = 1;
        zctype[i] = 'B';
        zindices[i] = i;
    }
    
TERMINATE:
    
    free_and_null ((char **) &protoObj);
    free_and_null ((char **) &protoSns);
    free_and_null ((char **) &protoRhs);
    free_and_null ((char **) &protoMatBeg);
    free_and_null ((char **) &protoMatCnt);
    free_and_null ((char **) &protoMatInd);
    free_and_null ((char **) &protoMatVal);
    
    if ( status ) {
        free_and_null ((char **) &zprobName);
        free_and_null ((char **) &zobj);
        free_and_null ((char **) &zrhs);
        free_and_null ((char **) &zsns);
        free_and_null ((char **) &zmatBeg);
        free_and_null ((char **) &zmatCnt);
        free_and_null ((char **) &zmatInd);
        free_and_null ((char **) &zmatVal);
        free_and_null ((char **) &zlb);
        free_and_null ((char **) &zub);
        free_and_null ((char **) &zctype);
        free_and_null ((char **) &zindices);
    }
    else {
        *numcols_p   = nCols;
        *numrows_p   = nRows;
        *numvars_p   = nVars;
        *objsen_p    = CPX_MAX;   /* The problem is maximization */
        
        *probname_p  = zprobName;
        *obj_p       = zobj;
        *rhs_p       = zrhs;
        *sense_p     = zsns;
        *matbeg_p    = zmatBeg;
        *matcnt_p    = zmatCnt;
        *matind_p    = zmatInd;
        *matval_p    = zmatVal;
        *lb_p        = zlb;
        *ub_p        = zub;
        *ctype_p     = zctype;
        *indices_p   = zindices;
    }
    return (status);
}


int createProtoType(double **protoObj_p, int **protoMatBeg_p, int **protoMatCnt_p, int **protoMatInd_p, 
                    double **protoMatVal_p, char **protoSns_p, double **protoRhs_p)
{
    
    int z,x,i,j;
    int nzConstI=2, nzConstII=1, row=0, status=0;
    
    int nzConsts         = nzConstI + nzConstII;
    int zCols            = nzConsts * nEdges;
    int nzVarsPerCluster = nzConsts * nEdges;
    int nxVarsPerCluster = nEdges * (nzConstI + nzConstII*2) + 2*nVertices;
    int nVarsPerCluster  = nzVarsPerCluster + nxVarsPerCluster;
    int nColsPerCluster  = nVertices+nEdges;
    int nRowsPerCluster  = nzVarsPerCluster + nVertices + 1;
    int totalNumOfConsts = nzConsts*nEdges + nVertices + 1;
   
    double *protoObj = NULL;
    char *protoSns   = NULL;
    double *protoRhs = NULL;
    int *protoMatBeg = NULL;
    int *protoMatCnt = NULL;
    int *protoMatInd = NULL;
    double *protoMatVal = NULL;
    
    int *X = NULL, *xCount = NULL, *xStart = NULL, *xConstsType = NULL;
    
    protoObj    = (double*)calloc(nColsPerCluster, sizeof(double));
    protoSns    = (char*)calloc(nRowsPerCluster, sizeof(char));
    protoRhs    = (double*)calloc(nRowsPerCluster, sizeof(double));
    protoMatBeg = (int*)calloc(nColsPerCluster, sizeof(int));
    protoMatCnt = (int*)calloc(nColsPerCluster, sizeof(int));
    protoMatInd = (int*)calloc(nVarsPerCluster, sizeof(int));
    protoMatVal = (double*)calloc(nVarsPerCluster, sizeof(double));
    
    X           = (int*)calloc(totalNumOfConsts, sizeof(int));
    xConstsType = (int*)calloc(totalNumOfConsts, sizeof(int));
    xCount      = (int*)calloc(nVertices, sizeof(int));
    xStart      = (int*)calloc(nVertices, sizeof(int));
    
    if (protoObj    == NULL || protoRhs    == NULL ||
        protoSns    == NULL ||
        protoMatCnt == NULL || protoMatBeg == NULL ||
        protoMatInd == NULL || protoMatVal == NULL ||
        X           == NULL || xConstsType == NULL ||
        xStart      == NULL || xCount      == NULL)  {
        
        status = 1;
        goto TERMINATE;
    }
    
    /* implementing Z's obj, matBeg, matCnt, matInd, matVal, sns, rhs */
    /* preparing X[], xConstsType[] for later usage to implement constraints I & II */
    for(z=0; z<nEdges; z++){
        
        protoObj[z] = getEdgeWeight(z);
        protoMatBeg[z] = row;
        
        for(i=0; i<nzConsts; i++){
            
            protoMatCnt[z]++;
            protoMatInd[row] = row;
            protoMatVal[row] = 1;
            
            X[row] = getVertexID(z, i);
            xConstsType[row] = (i == 2) ? 2 : 1;
            
            protoSns[row] = (i == 2) ? 'G' : 'L';
            protoRhs[row] = (i == 2) ? -1 : 0;
            
            row++;
        }
    }
    
    /* implementing X's obj, sns, rhs */
    /* preparing X[], xConstsType[] for later usage to implement constraints III & IV */
    for(x=0; x<=nVertices; x++){
        
        X[row] = x;
        xConstsType[row] = (x == nVertices) ? 4 : 3;
        
        protoSns[row] = (x < nVertices) ? 0-'E' : 0-'G';
        protoRhs[row] = 1;
        
        if (x < nVertices) protoObj[z+x] = 0;
        
        row++;
    }
    
    /* preparing xCount[] using X[] & xConstsType[] */
    for (i=0; i<totalNumOfConsts; i++){
        
        if (xConstsType[i] % 2){
            xCount[X[i]]++;
        }
        else{
            for(j=i-1; j>=i-X[i]; j--){
                
                xCount[X[j]]++; 
            }
        }
    } 
    
    /* implementing X's matCnt, matBeg using xCount[] */
    /* preparing xStart[] */
   
    for(x=0, i=0; x<nVertices; x++){
        
        xStart[x] = i;
        i += xCount[x];
        
        protoMatCnt[z+x] = xCount[x];
        protoMatBeg[z+x] = protoMatBeg[z-1+x] + protoMatCnt[z-1+x];
        
    }  
    
    
    /* implementing X's matInd, matVal using X[], xConstsType[], xStart[] */
    for (i=0; i<totalNumOfConsts; i++){
        
        switch (xConstsType[i]) {
            /* const I */
            case 1:
            {
                protoMatInd[zCols + xStart[X[i]]] = i;
                protoMatVal[zCols + xStart[X[i]]] = -1;
                xStart[X[i]]++;
                
            }
            break;
            /* const II */
            case 2:
            {
                for(j=i-1; j>=i-X[i]; j--){
                    
                    protoMatInd[zCols + xStart[X[j]]] = i;
                    protoMatVal[zCols + xStart[X[j]]] = -1;
                    xStart[X[j]]++;
                }
            }  
            break;
            /* const III */
            case 3:
            {
                protoMatInd[zCols + xStart[X[i]]] = -3;
                protoMatVal[zCols + xStart[X[i]]] = 1;
                xStart[X[i]]++;
            }
            break;    
            /* const IV */
            case 4:
            {
                for(j=i-1; j>=i-X[i]; j--){
                    
                    protoMatInd[zCols + xStart[X[j]]] = -4;
                    protoMatVal[zCols + xStart[X[j]]] = 1;
                    xStart[X[j]]++;
                }
            }
            break;
        } 
    }
    
    
TERMINATE:
    
    free_and_null ((char **) &X);
    free_and_null ((char **) &xConstsType);
    free_and_null ((char **) &xStart);
    free_and_null ((char **) &xCount);
    
    if ( status ) {
        free_and_null ((char **) &protoObj);
        free_and_null ((char **) &protoSns);
        free_and_null ((char **) &protoRhs);
        free_and_null ((char **) &protoMatCnt);
        free_and_null ((char **) &protoMatBeg);
        free_and_null ((char **) &protoMatInd);
        free_and_null ((char **) &protoMatVal);
        
    }
    else {
        *protoObj_p    = protoObj;
        *protoMatBeg_p = protoMatBeg;
        *protoMatCnt_p = protoMatCnt;
        *protoMatInd_p = protoMatInd;
        *protoMatVal_p = protoMatVal;
        *protoSns_p    = protoSns;
        *protoRhs_p    = protoRhs;
        
    }
    
    return (status);
}

